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The determination of the primary cosmic ray mass composition from the longitudinal development of atmospheric 
cascades is still an open problem. In this work we propose a new method of the multiparametric topological 
analysis and show that if both X max - the depth of shower maximum and N ma x - the number of charged particles 
in the shower maximum are used, reliable results can be obtained. 



1. Introduction 

The study of the longitudinal profile of indi- 
vidual atmospheric cascades started in the early 
eighties with the development of the fluorescent 
light detection technique implemented for the 
first time within the framework of the Fly's Eye 
experiment P^. After these pioneering efforts, 
there has been only one experimental array con- 
tinuing this type of studies: HiRes j2j and only 
recently a new and much more powerful detec- 
tor has started to collect data: the Fluorescent 
Detector (FD) of the Pierre Auger Observatory 
This instrument will produce a large data 
flow over the next decades and is therefore call- 
ing for new and accurate data analysis procedures 
capable to fully exploit the large amount of infor- 
mation contained in the FD data. 

It is rather surprising, that while there 
are many methods, both parametric and non- 
parametric (KNN, Bayesian methods, pattern 
recognition, neural nets etc.), used to discrimi- 
nate individual cascades on the basis of ground- 
based information, very little has been done to 
fully exploit the amount of information contained 
in FD data. To our knowledge, the only method 
available so far makes use of X max - the depth 
of the maximum cascade development 0] and de- 
rives the observed mean mass composition as a 
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function of the primary energy. Since there is 
a minimum bias in the detection of cascades of 
different origin, the observed mass composition 
coincides practically with the primary composi- 
tion. It has to be stressed also that this approach 
relies on statistical grounds and therefore does 
not allow the identification of the primary par- 
ticle for each individual cascade. Furthermore, 
even though in the longitudinal profile the X max 
parameter is the most sensitive to the mass of the 
primary particle, its sensitivity is still weak. For 
instance, at a primary energy of 1 EeV (10 18 eV) 
the mean iron induced cascade has X max only 
(11-12)% lower than for a proton induced one, 
i.e. a difference which is of the same order of 
magnitude as the intrinsic fluctuations of X max . 

As we shall discuss below, such unsatisfactory 
situation improves drastically if other, seemingly 
less significant parameters, are taken into ac- 
count. Among them, we might have: N max - 
the number of particles (mostly electrons) in the 
maximum of the cascade, the speed of rise in 
the particle number etc. For instance, at fixed 
primary energy, N max is about the same for all 
cascades, and even though iron induced cascades 
produce more muons and less energy is carried 
out by the electrons, the effects on N max are very 
small. However, due to the lower energy per con- 
stituent nucleon in the primary iron nucleus, the 
cascade development and the rise of the cascade 
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curve are on average faster than for cascades orig- 
inating from protons. This useful information is 
neglected when only N max is taken into account. 
This type of arguments triggered our efforts to 
find methods capable to fully exploit the infor- 
mation contained in the cascade curves and to 
allow the identification (at least in terms of prob- 
abilities) of the cascade origin also for individual 
showers. 

Although tailored for possible applications in 
the context of the Pierre Auger experiment, 
methods described below are quite general and 
may find application in other similar experiments. 
However, in what follows we shall focus on data 
similar to those expected from the Fluorescent 
Detector (FD) of the Pierre Auger Observatory: 
namely on the longitudinal profile of each atmo- 
spheric cascade, i.e. on the number of charged 
particle N c h as a function of the atmospheric 
depth X. As it will be shown below, this pro- 
file carries more information than X max or N max 
alone. 

We would like to stress that even though the 
fraction of the hybrid events, in which the infor- 
mation on the shower from the Pierre Auger Sur- 
face Array is supplemented by the Fluorescent 
Detector data, will hardly exceed 10% of the to- 
tal statistics accumulated by the Surface Array, 
these events need to be properly handled since 
they contain the maximum information. In this 
paper we restrict ourselves to the analysis of the 
longitudinal development of cascades. The treat- 
ment of hybrid data will be addressed in subse- 
quent papers. 

2. The simulated data 

In what follows, we assume that the primary 
energy estimates for hybrid events will be accu- 
rate at a few percent level The data set used 
to implement and test the methods described in 
the following sections consists of 7600 vertical cas- 
cades produced by particles with the fixed energy 
of 1 EeV, simulated using the CORSIKA pro- 
gram (version 6.015, 0) with the QGSJET in- 
teraction model. Simulations were performed at 
the Lyon Computer Centre. The primary nuclei 
were P, He, O and Fe, each of them initiating 



1900 cascades. The CORSIKA output provides 
the number of charged particles at atmospheric 
depths sampled with 5 gcm~ 2 intervals. 

We clipped the data at a depth of 200 g cm~ 2 , 
since the FD detection threshold does not allow 
to detect the weak signals at the beginning of 
the cascade development. The maximum atmo- 
spheric depth was set at 870 g cm -2 , roughly cor- 
responding to the level of the Pierre Auger Ob- 
servatory. In Figure 1 we show a subsample of 50 
cascades for each primary. 



Fi gure 1 . A few examples of the longitudinal develop- 
ment of 1 EeV vertical atmospheric cascades induced by 
protons P, helium He, oxygen O and iron Fe nuclei, spec- 
ified in headers. Each plot shows 50 typical cascades. 



3. Multiparametric Topological Analysis 
(MTA) 

The MTA method, when applied to FD data, 
relies on a topological analysis of correlations 
between the most significant parameters of the 
shower development in the atmosphere. In princi- 
ple this method could be used also with a greater 



x 10 

, 8000 

» 

j 7000 
[6000 

; 5000 

ft 

i 4000 

• 3000 
! 2000 
! 1000 
: 

x 10 

, 8000 

> 7000 

[6000 

J 5000 
n 

; 4000 

* 3000 
! 2000 
! 1000 

! o 



Prot 


mi, lcgE,eV= 18,0=0* 































200 400 600 800 

' X (g/cm 2 ) 



; Oxygen, log E,eV=1i 8,0=0? 































200 400 600 800 
X (g/cm 2 ) 



x 10 
, 8000 
> 7000 
[6000 

: sooo 

i" 4000 
1 3000 
! 2000 

! 1000 

! o 

x 10 
i 8000 
\ 7000 
[6000 
i 5000 
" 4000 
\ 3000 
\ 2000 
I 1000 
! 



Helium. logX,eV= 13.0=0* 






























3 



number of parameters, however, in this paper 
we restrict ourself to the simple case of two pa- 
rameters only: X max (the atmospheric depth of 
the shower maximum) and N ma2; (the number of 
charged particles at the depth of X max ). A scat- 
ter plot of these two parameters has been built 
using the 4 x 1000 showers already described. Fig- 
ure [3 shows the scatter plot for only proton and 
iron induced showers while their projected distri- 
butions are shown in Figure |3J 



Figure 2. N maa; vs. X ma:E scatter plot for proton (star) 
and iron (circle) induced showers 



The basic idea of the MTA method is to divide 
the entire area of the scatter plot into cells whose 
dimensions arc defined by an accuracy with which 
the parameters can be measured. In our simula- 
tions the value of 20 g/cm 2 was used as the width 
of X max bin, while 50xl0 5 was assumed for the 
width of ~N max bin. In each cell we can find the 
total number of showers N| ot , as the sum of N P , 
Njy e , No and N Fe showers induced by P, He, O 
and Fe respectively, and then derive the associ- 
ated frequencies: p P =N p /Nj ot , p}j e =W He /W tot , 



Proton (block/dashed) — Iron (blue/full) 

(O 400 F 




Figure 3. Projected distributions of the scatter plot in 
Figure [5] for protons (dashed line) and iron nuclei (full 
line) 



Po =N /Nj ot and pj, e =Nk/Nj ot which can be 
interpreted as the probability for a real shower 
falling into the i th cell to be initiated by pro- 
ton, helium, oxygen or iron primary nuclei. In 
other words, in the case of an experimental data 
set of N exp showers, it may be seen as composed 
of a mixture of N exp x pp proton showers, N exp 
x pue helium showers, N exp x po oxygen show- 
ers and N exp x pp e iron induced showers, where 
PA = £ 4 pV 

In order to generate the X max — N max scatter 
plot and produce the relevant matrix of cells we 
used a set of 4x1000 simulated showers. Then 
we used another subset of 4x300 showers to de- 
termine the probabilities Py . For each individual 
shower i in a given subset the partial probabili- 
ties p P , p\j e , Po and p l Fe have been read from 
the relevant cell i. The sum of such probabili- 
ties over the entire subset of 300 showers permits 
to estimate the probability pa for a shower of a 
given nature A to be identified as a shower gen- 
erated by the P, He, O or Fe primary particle. 
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This probability is shown in Figure One can 
see that the method works quite well being able 
to attribute the highest probability to the correct 
nuclei. 



Figure 4. Application of the MTA method to the N max - 
y<-max scatter plot - mean probability for the cascades 
induced by i nucleus: P, He, O and Fe (specified in head- 
ers) being identified as induced by j nucleus, indicated at 
abscissa . 



4. Determination of the primary mass 
composition 

The obtained mean probabilities Py for cas- 
cades induced by type i nuclei to be identified as 
those induced by j nuclei for the pure primary 
mass composition can be used for the reconstruc- 
tion of the mixed primary mass composition as 
the coefficients in the system of linear equations: 

n' P = Y,1 =l n Ai ■ PAi-,p 
n-He = ^f =1 n Az ■ P A ^He (1) 

n O = ^i=l n A t ■ PAi^O 



n Fe = ^i=l n Ai ■ PAi^Fe 

where n Ai (np, riHe, no and nFe) are the true 
numbers, which determine the primary mass com- 
position in the sample of N = np+nne+no+np e 
cascades, which are observed as n' P , n' He , n' and 
n' Fe , due to a misclassification. 

In order to invert the problem and to recon- 
struct the abundances n Ai in the primary mass 
composition from the observed abundances n' A . 
and the known probabilities Py, we can apply 
any method capable to solve the inverse problem 
taking into account possible errors of the observed 
distribution and the constraint: 

ZUn Ai = = N (2) 

The observed abundances were simulated using a 
set of 4x600 cascades different from those used 
for the determination of Py . At the moment 
the result of the solution for the uniform primary 
mass composition with all abundances p A = 
= 0.25 reconstructed by MTA method is pp = 
0.25, pHe = 0.24, po = 0.23 and p Fe = 0.28. 

As an example of the non-uniform mass com- 
position we have taken the case with pp = pp e 
= 0.1, pHe — Po = 0.4, which resembles the 
mass composition of cosmic rays at the knee en- 
ergy of ^3 PeV 0. The total number of cas- 
cades used for the simulation of non-uniform mass 
composition is 1500 (150P + GOOHe + 6000 
+ 150^). The result of the solution for this 
case is pp = 0.09, pu e = 0.40, po = 0.35 and 
PFe = 0.16. 

It is seen that absolute values of the abundance 
are reconstructed satisfactorily. For the time be- 
ing the values of the probabilities are assumed 
to be known precisely. A more accurate solution 
of the system, taking propely into account the 
errors of the probability, is under study. 

5. Application for experimental cascades 

It has to be stressed that all the results outlined 
above are biased by the fact that we are dealing 
with simulated 'noiseless' data and more realistic 
testing has to be performed using data which take 
into account the varying primary energy, the incli- 
nation angle, the instrumental signature (noise) 
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and the various sources of errors (uneven and in- 
complete sampling, etc.). Application of MTA is 
straightforward - the relevant simulations should 
be made with a maximum possible statistics. 

6. Conclusions 

We proposed and tested the Multiparametric 
Topological Analysis method for the determina- 
tion of the mean primary cosmic ray mass compo- 
sition on the basis of measurements of the longitu- 
dinal development of atmospheric cascades. The 
method employs more information about the lon- 
gitudinal development of atmospheric cascades 
than just the depth of the cascade maximum. 
Definitely the method needs further developments 
before being used for processing real experimental 
cascades, but the first results are very encourag- 
ing. The next efforts should include the increase 
of simulation statistics, the study of the effect of 
the shower energy spectrum and angular distri- 
bution, the accurate estimate of the reconstruc- 
tion errors. The great advantage of the proposed 
MTA method is that it is extremely easy to use 
and generalize for a larger number of observables. 
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